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1. Introduction 

The zero-temperature random-field Ising model (zt-RFIM) is a prototype model for 
systems that exhibit avalanche dynamics when slowly driven accross athermal first- 
order phase transitions [H [2] . Examples of such transitions include the condensation of 
fluids in porous media [3], the martensitic transformation or magnetisation reversal 
of ferromagnets (Barkhausen effect) [5]. For experimentally reasonable time scales, 
thermal fluctuations do not play an important role and the dynamics in this kind 
of system proceed along a non-equilibrium path consisting of metastable states [6]. 
Avalanches are a manifestation of such behaviour and correspond to the driven-induced 
passage between two metastable states. The zt-RFIM predicts that the properties of 
avalanches (e.g. their size or duration) are drastically affected by the degree of quenched 
disorder in the system. More explicitly, the model predicts three different regimes for 
avalanche behaviour depending on the degree of disorder. The main feature of the small- 
disorder regime is the existence of an avalanche that is infinite in extent. In contrast, 
high degrees of disorder lead to a regime where all avalanches are small. These two 
regimes are separated by an intermediate situation where the model exhibits critical 
behaviour (i.e. a continuous phase transition) [H El El E] . For any degree of disorder, 
avalanche- like dynamics lead to spacial correlations in observable properties averaged 
over quenched disorder (e.g. magnetisation, stress or fluid density), which extend over a 
typical length scale, called the correlation length. This quantity is expected to diverge 
in the critical regime, and this has been shown to be the case in the systems described 
by mean field [7] and hypercubic lattices [9]. Concerning other topologies, including a 
Bethe lattice, the situation is less clear. 

Previous studies of the zt-RFIM on the Bethe lattice have obtained exact results 
for the magnetisation [8] , the avalanche size distribution [10] , the different contributions 
to the energy [TT], and the number of metastable states [121 lEl HI]- In particular, 
it has been demonstrated that for coordination numbers g > 3 the system exhibits 
a discontinuity in its magnetisation hysteresis loop for small amounts of disorder 
associated with the infinite avalanche. The universality class of the critical point on 
such a lattice has been suggested to be the same as that of the mean-field system 
(complete graph) [H [7], with the critical exponents for the order parameter, /3 = 1/2 
and 5 = 3, being identical to the mean- field values [H [15]. The spatial correlations 
of the model have also been investigated in the past but mostly using numerical 
simulations [9] or approximate analytical methods based on mean-field descriptions 
and/or renormalization group (RG) techniques [16]. In this paper, we derive an exact 
analytical expression for the spin-spin correlation function corresponding to the zt-RFIM 
with spins placed on a Bethe lattice. Our results are exact for any coordination number 
and confirm the validity of the functional form for the correlation function derived 
independently in [17], where it is assumed that the extension from the one-dimensional 
case (i.e. with g = 2) is valid. 

It is well established that in a hypercubic lattice of dimensionality d around 
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criticality, the correlation function, C{r), decays exponentially with distance r for the 
zt-RFIM, 

C(r) = A(r)exp(-r/0 , (1) 

where A{r) obeys a power law for large r, and ^ is the correlation length. At the critical 
point, ^ — )■ oo, and so A{r) represents the critical behaviour of the correlation function, 
found to be 

A{r) oc l/r'^~^+^ , 

with anomalous dimension fj [71 [H]. In a prototype loopless topology, the Cayley tree 
or Bethe lattice, the correlation function exhibits additional exponential behaviour due 
to the associated hyperbolic topology, 

Cir) = B{r)iq-l)-^expi-r/0 , (2) 

where q is the coordination number and B{r) is a power-law function for large values 
of chemical distance r. The exponential prefactor (g — 1)'' in ([2]) gives the number of 
lattice sites within r shells of the Bethe lattice, and plays the role of the factor r"^ in the 
equivalent expression for hypercubic lattices. The function 

B{r) oc , 

with some exponent fj (found below to he rj = —1), accounts for all power-law 
behaviour of the correlation function. There is some inconsistency in the definition of 
the correlation length for two-state models (the equilibrium Ising model, the zt-RFIM 
and percolation) on the Bethe lattice. In the equilibrium Ising model [18] and the zt- 
RFIM [17] on the Bethe lattice, the prefactor (g — 1)~'' in ([2]) has been absorbed into 
the definition of the correlation length 

C(r) oc (g — 1)^^' exp(— r/^) = exp(— r/ro) , (3) 

with 

ro = (ln(g-l)+r')"' , (4) 

chosen to be the correlation length. In an alternative definition, for both percolation [19] 
and the Ising model [20], the correlation length has instead been identified with ^ in (j2]), 
thus separating the length scale of dynamic correlations from the intrinsic topological 
length scale of the Bethe lattice, l/ln(g — 1). 

Below, we demonstrate that (i) the value of ^ diverges at the critical point according 
to a power law, in contrast to ro which remains finite at criticality and (ii) the critical 
exponents governing the divergence of ^ around the critical point are consistent with the 
mean-field values for the divergence of correlation length in hypercubic lattices above 
the upper critical dimension, d > dc = Q. The exponents, describing the divergence 
of C, in terms of the external field H and degree of disorder A, ^ oc {H — Hc)~'^ and 
^ oc (A — Ac)~'^, are found to be /i = 2/3 and z/ = 1. As argued below, these values can 
be related to those in the mean- field model, fi/d = /imf = 1/3 and v/d = pmy = 1/2, 
where d = 2 originates from fractal dimension of a random walk [191 El]- Therefore 
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these findings motivate us to choose ^ as the correlation length describing the critical 
behaviour of the zt-RFIM on a Bethe lattice. Such a definition, while in contrast with 
that used in [17], separates the topological contribution of the Bethe lattice to the 
correlation length, and allows its comparison with the correlation length for hypercubic 
lattices. 

2. Model 

The zt-RFIM involves a set of N Ising spins, {si = ±l|i = 1,2,...,A^}, interacting 
ferromagnetically with strength J > 0, affected by an external field H and quenched 
local disorder, manifested in independent random local fields, hi, at each site i, according 
to the following Hamiltonian, 

'H = -J^SiSj - Hj2si-J2^i^i ^ (5) 

(ij) i i 

where {ij) denotes a sum over all pairs of nearest neighbours. The random fields hi 
are assumed to be identically distributed according to the probability density function 
p{hi) with, for convenience, zero mean and standard deviation A. It has been found 
that in certain topologies of the network of spins, and for disorder A < A^, the zt- 
RFIM exhibits a discontinuity in the hysteresis loop for magnetisation, m = N^^ Yl,i Si, 
at some coersive field -ff(A), where the magnetisation in a large fraction of the system 
reverses in a single avalanche. Increasing disorder is found to reduce the size of this 
discontinuity, and to remove it continuously at the critical value of disorder Ac, around 
which spin-spin correlations become infinite-ranged and the system exhibits scale-free 
universal properties [H El [9], |22] • 

We start with a description of the dynamical rules governing the relaxation of spins 
within the zt-RFIM. The external field is initially set to H = — oo, forcing all spins to be 
Si = —1, and then allowed to increase adiabatically to = oo. When the external field 
is varied, the system becomes unstable and relaxes into a new metastable state through 
a series of spin flips. We assume that relaxation takes place according to Glauber single- 
spin flip dynamics for the zero temperature case [8], such that a spin flips only if such 
a flip reduces the overall energy. If a spin at node i changes state (i.e. Si — )■ — Sj), it 
induces a change in energy 2/jSj, where fi = JJ2{i\j) Sj + H + hi is a local field. If /jSj < 
then the spin flip will occur, otherwise it remains aligned with the local field. As spins 
flip one at a time (HE], the local field at surrounding sites changes, and other spins may 
become unstable and also flip, therefore relaxation occurs in an avalanche like manner. 
At a given external field, this process continues until all spins are aligned with their 
respective local fields at which point the system becomes stable, thus creating clusters 
of flipped spins around each nucleation site. As the external field increases avalanches 
continue to progress, and it is known that the set of spins that flip in multiple avalanches 
when the external field is slowly increased from — oo to if coincides with the set of spins 
that would flip if the field had stepped from — oo to H (see [8]). Accordingly, it can be 



Exact spin-spin correlation function for the zero-temperature random-field Ising model5 






h-2 

h-1 

g = h 



Figure 1. Caylcy tree of h generations, labelled from g — 1 to g = h, with generation 
g = 1 being the central chain (in bold at the top of the tree) of r sites used for 
calculation of the correlation function. In the infinite limit, the deep interior of this 
form of the Cayley tree has the same properties as the more standard form, where 
there is a single central (root) site, rather than a chain. 



assumed that the correlation function calculated at any external field depends only on 
the properties of avalanches occurring at that field. 

3. Spin-spin Correlation Function 

The spin-spin correlation function, C{r), between two spins labelled i = and i = r, 
separated by the shortest chemical distance r > 0, is given by, 

C{r) = {SoSr) - (So) {Sr) , (6) 

where (...) denotes an average over quenched disorder. In order to evaluate the 
correlation function, it is convenient to define a Cayley tree of height h by dividing 
it into a set of h generations g of spins (see figure [1]). The first generation, g = 1, forms 
a central chain consisting of r + 1 connected spins, in contrast to the standard definition 



with a single central (root) site (presented in Appendix A). The two boundary spins of 



the central chain, i = and i = r, are each connected to one spin within the central 
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chain, and to g — 1 spins in generation g = 2 oi the Cayley tree. The interior spins of the 
central chain < i < r interact with 2 other neighbouring spins (see thick horizontal 
lines in figure [I]) in the central chain and q — 2 second generation spins. Each spin in 
generation g > 1 then interacts with a single spin in generation — 1, and q — 1 spins in 
generation g + 1. Spins in the last (boundary) generation, g = h, interact only with a 
single spin in generation h — 1. In the limit of large number of generations, a Cayley tree 
defined in this way tends to the same limit, i.e. a Bethe lattice, as the more standard 
Cayley tree defined in [HI ESj . 

For finite values of external magnetic field, spins in all generations are relaxed. 
First, we relax the spins in generations g > 1, so that the spins in generation g = 2 are 
up with probability P*, which is the solution to the self-consistent equation, 

P* = F{P*) . (7) 

The function F{P*) is given by [8], 

li-H + (q- 2m) J, oo) , (8) 

with, 

h" 

I{h',h")= Jp{hi)dh. (9) 

h' 

Then, knowing the state of spins in generation g = 2, we relax the spins in the central 
chain g = 1 and calculate the correlation function (see Appendix A). 

The relaxed central chain can contain both fiipped and unfiipped spins. In order 
to know the state of spin i G [0, r] we need to know the local field, fi, which depends on 
the neighbourhood of this spin and the random field, hi, 

/, {h„ n[, O = H + h,- (2(n: + <) - q)J . (10) 

Here, the state of the neighbourhood is represented by the number n[ of fiipped 
neighbours in the central chain and the number n'- of fiipped neighbours in generation 2 
prior to the relaxation of the central chain, so that the total number of fiipped neighbours 
= n- + n'-. 

According to the values of the local fields, fi {hi,n[,n'-), where the variables hi 
and n'- are fixed during the relaxation of the central chain, spins in nodes < i < r 
can be divided into three categories: (1) those that experience a positive local field 
when all neighbours within the central chain are down n'^ = 0, i.e. spins in the 
set Ri = {{n",hi)\fi{hi,0,n") > 0}; (2) those that only experience a positive 
local field after one of their two neighbours in the central chain has fiipped i.e. 

= {{n'l,hi)\{fi{hi,0,n'{) < and f; (h;, 1, nf) > 0)}; (3) all other spins R3 = 
{«,/i.)|/.(/i„l,<)<0}. 

In general, the states of the chain-boundary spins i = and i = r are determined 
entirely by the categories q of all the spins in the central chain. Therefore, one can 
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calculate the value of (sos^) given the probabilities Qc = Qdl — 2) and Q'^ = Qdq — 1) 
of interior and boundary spins in the central chain falling into category c = 1,2, or 3, 
where, 



-'5 



QM 




mill) ^max) j (H) 



and /imin and /imax 

q], oo), c=l 

(^min, /^max) = <! J[2n"-q]), C=2 ^^^^ 

J[2n" + 2 - g]) , c=3 

The value of Qdq) given by (fTTj) is a sum of probabilities P{n") that n" neighbours 
in generation g = 2 have flipped, multiplied by the conditional probability, 
P {{n", hi) G Rc\n"), that the spin is in category c given this value of n". This conditional 
probability translates into the integral /(/imin, ^max) over the range of random fields 
specified by (fT2l) which cause the spin to be in that category. 

Avalanches within the central chain are nucleated at all sites of category 1, and 
propagate through sites of category 2, terminating either at the end of the chain, or at a 
site of category 3. Therefore, for a given configuration of categories Cj of spins in nodes 
< i < r it is possible to determine the state of the spins at the end of the chain, i = 
and i = r, by finding the first site from each end which is not of category 2. If this first 
site not of category 2 is of category 1 then the respective end spin will be up and if such 
first site is of category 3 the end spin will be down (see Appendix B). Summing over all 
combinations of such systems it is possible to calculate {soSr). Similarly, it is possible 
to consider just one spin (sq). Therefore the correlation function can be found by the 
formula IQ, and is given by. 



C{r) 



(i-Q2m 

2Q'2- 



Q'2 



+ 7'^dl.2 r ) . (13) 



'i-Q2)Ql' 



\ Q2, 

Equation f[T51) is the main result of our derivation, and represents the exact expression 
for the correlation function in a Bethe lattice. For r ^ 1, the correlation function should 
follow the behaviour [201 



CW~7 ^expl — 1 . (14) 




(g - l)'"r'? 

where, by comparing (0311 with ( lT3l) in the limit of large r, it is found that fj = —1 and 
the correlation length is given by, 

e = [-ln((g-l)g2)]"' . (15) 
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Figure 2. The dependence of the correlation function C(r) on r for g = 4 
according to (lines) and numerical simulations (symbols) for several values of 
external field: H = 0.75 (circles, solid line), 1.0 (squares, dashed line), and 1.25 
(diamonds, dot-dashed line). All results were obtained for systems characterised by 
A — Ac = 1.78125895. Numerical simulations consider a g-regular graph {q = 4) of 
N — 10^ spins, and calculate the mean of the correlation function between all pairs of 
sites in the system for 10"^ realisations of quenched disorder. 



The exact expression for correlation function given by ( fT3l) can be supported 
numerically by evaluation of the correlation function on random g-regular graphs. These 
are random graphs with a fixed coordination number q whose topology at local scales is 
similar to that of a Cayley tree with the same coordination number [23]. The absence 
of boundaries in g-regular graphs make them more suitable for numerical simulations 
than the Bethe lattice used in our analytical derivation (see Appendix C). Figure [2] 
demonstrates that the results of numerical calculations closely match the analytical 
formula, (|T3l) . 



4. Expansion Around Criticality 



Near the critical point the correlation length diverges according to standard scaling 
relations, 



e oc (A - Ae; 



if H : 
if A 



He 
^ A. 



(16) 
(17) 



with V and ^ being the critical exponents. In order to find these exponents, it is 
convenient to rewrite (fTS!) using the following property of Q2-, 

1 dF{P*) 



(g-1) dP* ■ 
Introducing this expression in f[T^ . one obtains, 

dG 



In 1 + 



dP* 



:i8) 



(19) 
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where G{P*) = F{P*) — P* . Equation ( TT^ can be expanded in terms of dG/dP* which 
is equal to at criticahty and small nearby, resulting in, 



In the case of Bethe lattices with coordination numbers g = 2 or g = 3, there is no 
transition at any non-zero disorder. However, for such coordination numbers, it is 
possible to show by our method that the correlation length diverges exponentially as 
disorder goes to zero. This result is in agreement with the calculation given in [25] for 
the case q = 2. 

For higher coordination numbers g > 4, where a transition exists, expansion of (!20ll 
around criticality can be done in terms of H — He, A — and P* — P*. Using the 
known mean- field exponents S = 3 and /3 = 1/2 [H E] describing the behaviour of 
P* — Pc in terms of H — He and A — Ac, respectively, around criticality and noting that 
both the first and second derivatives dC/dP* = d^G/dP*^ = one obtains the values 
= i^MFd = 1 and fi = fiupd = 2/3 for the correlation length critical exponents (see 
Appendix D). 

The values of these exponents can be explained in terms of the mean-field values 
using the argument for percolation of [191 [21]. this argument, it is demonstrated that 
the critical percolation cluster in high dimensions {d > 6) is sparse enough as for the 
loops to be unimportant. In this case, the critical percolation cluster is essentially the 
same as that on a Bethe lattice, i.e. its backbone consists of chains of links (which form 
random walks), with occasional branch points (nodes). The mean chain length between 
branching nodes is defined as the correlation length which scales as {p — Pc)~'^'^^"' (where 
p and Pc are bond probability and the critical bond probability) with the power Uperc = 1 
for chemical distance and Uperc = 1/2 for Euclidean distance. The same arguments can 
be applied to the critical avalanche clusters in the zt-RFIM on a hypercubic lattice and 
Bethe lattice, so that the correlation length should scale with the exponents = 2/3, 
V = 1 for chemical distance through the avalanche cluster and /z = 1/3, v = 1/2 for 
Euclidean distance across the hypercubic lattice. 

5. Conclusions 

To conclude, we derived an exact analytical formula for the spin-spin correlation 
function for the zt-RFIM defined on a Bethe lattice and investigated the scaling of 
the correlation length near criticality. We demonstrate that if the correlation length 
is defined not to include the intrinsic topological length scale, then it diverges around 
criticality with exponents consistent with the mean-field description. In fact, choosing 
the dimensionality of the Bethe lattice to be a function of the length scale r [201 123] , 
that is, replacing the term r'^ in the standard definition of the correlation length with 
{q — ly, allows for analysis consistent with that performed for hypercubic lattices. 




(20) 
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i = i = 1 z = r — 1 i = r 




Figure Al. Caylcy tree of height h, drawn with each row of spins representing a 
generation g, 1 < g < h. The first generation, being the top of the tree, is a central 
chain of length r. The arrows in generations g = h, g = h— 1 and g = h — 2 represent 
the state of the spins after the second relaxation step, i = 2. In generation g = h, 
spins marked with | remain in the initial down-state after the £ = 2 relaxation, spins 
marked with f flipped on the first step of relaxation, £ = 1, and spins marked with 
tl" flipped on the second step of relaxation, £ = 2, owing to the increased local field 
caused by the spin flips occurring directly above them. 
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Appendix A. Relaxation of Generations g > 1 oi the Cayley tree 

In this section, we describe the relaxation of spins in all generations with g > 1 using 
a procedure similar to that described by [8] for obtaining the magnetisation in the 
zero-temperature random field Ising model (zt-RFIM). The system of spins within 
the ferromagnetic zt-RFIM, relaxing according to zero-temperature Glauber dynamics, 
obeys the Abelian property [1], i.e. the final metastable state of the spin system after 
relaxation is independent of the order in which spins relax. We therefore choose a 
relaxation order convenient to analysis. We start relaxation with boundary generation 
g = h dX the bottom of the tree (see figure lAll) . At the first step, the spins in all 
generations except the boundary generation h are held in the their initial down (J,, 
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s = — 1) state, while allowing any spin in generation h to flip (change state, s — > —s) 
if the local field acting on that spin is positive, and thus the flip reduces the overall 
energy (see the up spins, in figure lATl in generation g = h). At the next step, spins 
in generation h — 1 are similarly allowed to flip if doing so reduces the overall energy 
(see the up spins, in figure \K\\ in generation g = h — 1), while the sites in generations 
g < h — 1 are again kept fixed down. Spins flipping in generation h — 1 will cause a 
change in the local field acting at their non-fiipped neighbours in generation h, which 
may become positive and thus cause these spins to flip as well (see the up spins, fl', in 
figure IXn in generation g = h). Such secondary flips in generation h cannot influence the 
spin in the above generation, which is already flipped and thus in its final state (at zero 
temperature). Because of the loopless structure of the Cayley tree, these secondary flips 
are isolated from the above generations. Therefore, they cannot affect those generations 
and, in particular, the central chain (^f = 1) in which we are interested for calculating 
the correlation function. At subsequent steps, relaxation progresses up the generations 
in a similar way until the spins in generation g = 2 are relaxed. 

The above relaxation procedure is a random process, due to the contribution of 
quenched independent random fields {h\^^\g = l,...,h} to the local fields {fi^^\g = 
1, . . . ,h} at generation g. At the first step, all other contributions to the local fields 
{fi^^} are identical for all sites, meaning that the local fields are independent identically 
distributed random variables. As a consequence any spin Sj in generation h may flip 
independently of all others in this generation with probability P^^\ the probability that 
//'^^ > 0. At the second step, randomness in the local fields at sites in generation h — 1 
is introduced by the random fields {/i^^"^^} at those sites in combination with the field 
produced by a random number of flipped neighbours in generation h, both of which are 
independent random variables. Therefore, all the flips in generation h — lai the second 
step are independent events, occurring with probability P^^^^\ Similarly, according 
to this relaxation procedure, at the i-th step, the flips in generation h — i -\- 1 are all 
independent events. 

In order to calculate the values of the probabilities P^^\ consider flrst an arbitrary 
spin Si in generation g < h and calculate the conditional probability Pn of it flipping 
given that it is surrounded by n spins in the up state. The condition for the spin to flip 
is that the local fleld, 

f^^^ = H + h^^ + jJ2s,, (A.l) 

m 

is greater than zero, i.e. 

f^<^\h\'^;n,) = H + h'f^ - J{q - 2n^^) > . (A.2) 

Here, rij-^^ = Ui = + Sj)/2 gives the number of neighbouring spins j (in 

generations ± 1) of the spin i (in generation g) which are in the up state, Sj = 1. 
The probability of the random fleld h^ meeting condition flA.2l) is, 

/oo 
e {fiihu n)) p{hi)dhi = I{-H + (g - 2n)J, oo) , (A.3) 
-oo 
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where Q{x) stands for the Heaviside step function, and I{h^i^, h^g^y) is given by 

The variable n in f lA.SP is a random variable depending on the orientation of 
neighbouring spins, which are located in generations g± = ± 1. For 1 < g < h, when a 
spin in generation g is first relaxed (at relaxation step h — g the neighbouring spin 
in generation g^ is in the down-state, while the q — 1 neighbouring spins in generation 
g^ may each independently (as stated above) be in the up-state with probability P^9+^\ 
Therefore, the total number of neighbouring spins which are in the up-state is binomially 
distributed according to n ~ B{q - 1,P(3+1))^ i.e. the probability of n neighbours of a 
spin at generation g being up is given by 

1' 



F{n,g) 



n 



-1-n 



The probability P^^^ can then be found as a sum over all possible configurations of n 
neighbouring spins in the up-state (0 < n < g — 1), 

. n 



q-l 



pia) = Y^F{n,g)p^ = Y. 



1 _ p(a+l) 



Pn 



(A.4) 



n=0 n=0 

which is a recursion relation for P*^^) valid for 1 < g < h. Far from the boundary 
of very large Cayley trees {1 < g <^ h, h ^ I), P^^^ tends to some limiting value, 
P* = limh^oo P^^\ identical for all interior generations. This leads to the self-consistent 
equation for P* (see ([7]) and ([8]) 



F{P*] 
p* 




'[1-P 



*] q—l—m 



Pv 



(A.5) 
(A.6) 



so that P*^^-* = P* for all g > 1, i.e. after the relaxation of all spins except those in the 
central chain, spins in generation 2 are in the up-state with probability P*. 

In the case of the standard Cayley tree (see figure lA2|) . the relaxation procedure 
described above allows the magnetisation of the central spin (generation g = 1 in 
figure IA2p to be calculated in the following way. As mentioned above, f lA.5|) and (]A.6P 
are only valid for 1 < g < h, and the central spin at generation g = 1 should be 
treated separately. The central spin has q neighbours in generation g = 2, which, on 
the relaxation of that spin, are in the up-state with probability P'-^^ = P*. The number 
of neighbours of the central site that are in the up-state is binomially distributed, 
n ~ B{q,P*) and the expression for P*^^) is given by (]A.4p with q — 1 replaced by q 
and P(9+i) by P*, i.e. 

q 



p(i) 




P 



Pn ■ 



(A.7) 



n=0 



The value of P^^^ calculated by solving the self-consistent equations (lA.SP and (]A.6P and 
substituting into (lA.Tp allows the mean magnetisation of the lattice to be evaluated, 
(m) = 2P'^^^ — 1. The result of this calculation for magnetisation vs external field is 
given in figure R3| showing the known spinodal transition at low disorder [8]. 
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Figure A2. Diagram of a Cayley tree with a single node (r = 0) in the first generation, 
5 = 1, as considered by jH]. 




Figure A3. Mean magnetisation, (m), vs external field, H, scaled by interaction 
strength, J, for zt-RFIM on a Bethe lattice of coordination number q = 4, as 
magnetic field is swept upwards from — oo to +oo. The solid, dot-dashed and 
dashed lines correspond to different degrees of disorder around critically: A = 1.65, 
A = Ac = 1.78215895 and A = 1.9, respectively. The critical field {H^ = 1.0) is 
marked by the dotted vertical line. In the region where the magnetisation curve is 
multi- valued (marked by circles), the system follows the lower branch, so that a jump 
in magnetisation (represented by the vertical arrow) occurs when the external field 
passes the coercive field strength. 
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Co = 22 2 2a = /32 2 2 2 2 c^ = 2 



fc=4 Z=6 
Figure Bl. Distinct configurations of spin categories in tlie central chain, as used 

in (Ell). 



Appendix B. Derivation of Correlation Function From Categories 

In order to calculate the correlation function, we need to know the states of the chain- 
boundary spins i = and i = r which are determined entirely by the categories q of all 
the spins in the central chain. However if the categories of the boundary spins are either 
1 or 3 (see section [3] for definition), their state depends only on their own category. In 
particular, if the category of a boundary spin is 1 then it will have a positive local field 
and fiip. Conversely, if a boundary spin has a category 3 then it will have a negative 
local field and will never fiip. In the case that a boundary spin is of category 2, its 
state depends on the categories of the other spins in the central chain. Specifically, a 
boundary spin (i = or i = r) of category 2 will fiip only if its neighbour in the chain 
(either spin z = lorz = r — 1) fiips first. The state of the spins i = 1 and i = r — 1, 
given that their respective neighbours at the boundary are of category 2, can be found 
in the same way. Therefore, the state of spin i = depends on and coincides with the 
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state of the first spin, A; (0 < /c < r), counting from the i = boundary, that is not of 
category 2. The state of spin k is, itself, determined only by its category a = c^, with 
Cfc = 1 or Cfc = 3 but, by definition, a 2. Similarly, the state of spin r is determined by 
the state of spin r — I {0 < I < r — k), the first spin counting from the i = r boundary 
which is not of category 2, being instead of category /3 = Cr-i. In the special case of all 
the spins in the central chain being of category 2, the values of a, P, k and I cannot be 
defined. However, this situation is simpler because all spins in the central chain after 
relaxation are necessarily in the down-state, Si = —1. 

The probability, P{k, a; I, (3), that all spins i {0 < i < k) are of category 2 with spin 
k being of category a, and simultaneously all spins j {r — I < j < r) are of category 2 
with spin r — l being of category (3, can be calculated explicitly (see figure IBTI for details). 
In figure IBll we show the central chain of length r with white circles representing spin 
of either category 1 or 3, and grey circles are used for spins of category 2. Circles 
with crosses represent spins for which the category is irrelevant to the state of the 
boundary spins. All distinct configurations of spin categories (excluding the special 
case when all spins are of category 2) are represented by rows (a)-(g). In configuration 
(a), the boundary spin are not of category 2, so that the categories of all other spins are 
irrelevant. In configurations (b) and (d), one of the boundary spins is of category 2, and 
the first spin which is not of category 2 linked to this spin does not coincide with the 
other boundary spin. In configurations (c) and (e), all the spins are of category 2, except 
one of the boundary spins, meaning a = (3. In configurations (f)/(g), both boundary 
spins are of category 2, and the first spins not of category 2, counting from spin i = 
and from i = r do not coincide/do coincide. There are two important observations to be 
taken into account in order to evaluate P{k, a; The first observation is that a and 
(3 are independent variables only if A; + Z < r (see configurations (a), (b), (d) and (f) in 
figure IBT|) . while otherwise, i.e. when k-\-l = r, they are equal, a = (3 (see configurations 
(c), (e) and (g)). The second point to consider is that the configurations where one or 
both of the boundary spins are not of category 2 should be treated separately (i.e. in 
configurations (a), (b), (c), (d) and (e)), because in such configurations the probability 
distribution of a and /3 are described by Q'^ and Q'^ rather than Q„ and Q/?- Bearing 
these two points in mind, the probability P{k, a; I, /3) for each of the configurations 
presented in figure IBll is given by the following expressions, 

Q'aQ'p k = l = (a) 

Q'aQ2Q2~'Qp k = 0, 0<l<r (6) 

Sa,pQ'aQ2Q2-' k = 0, I = r (c) 

P{k,a;l,/3) = { Q'pQ'2Q2~^Qa I = 0, < k < r (d) (B.l) 

5a,pQ'pQ'2Q2-' 1 = 0, k = r (e) 

Q'^'^Qi^'-^QaQp k>0, l>0, k + l<r (/) 
5a,f3Q2^Q'2^^"^Qa k > 0, I > 0, k + I = r (g) 
The probability of each configuration is given by a product of probabilities that each 
relevant spin is of a specific category. For example, for configuration (g) the quantity 
refers to the probability of both boundary spins being of category 2, the power 
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Q 



l+k-2 



refers to the probability oil + k — 2 interior spins also being of category 2 and 



the term describes the single interior spin not of category 2. The probability of the 
special case that all spins are of category 2 is Q'2Q2~^- 

The states of spins and r are determined by a and (3 in (IB.ip . so that Sq = S{a) 
and Sr = S{f3) where, 



S{a) 



1, a = 1 
—1, a = 3 



Then the expectation value can be calculated by the sum, 

(sosr) = Y.i2'T.S{a)S{/3)P{k,a-l,(3) + Q','Ql-' 

a,/3 k=0 1=0 

Substitution of the probabilities given by (IB.ip into (1B.3P gives, 
(sosr) = {Q[ - Q'sY + 2{Q[ - Q',)Q', J2 QV\Qi - 



(B.2) 



(B.3) 



i=l 



■r-1 



+ g'/ ^(z - i)Qr'(Qi - Q^f + Q'2\r - m,'\i - Q2) 

i=2 

+ 2g^gr'(i - Q'2) + Q'2Ql''Q'2 ■ (b.4) 

The value of (sq) = {sr) in (I6l) can be found similarly by considering a central chain 
of infinite length (r = 00) and evaluating the probabilities, P{k,a), that the sites i, 
< i < k {k > 0) are of category 2 and that spin k is of category a, 

^^^'"^^{ g^gru, k>o ■ 

Therefore the mean magnetisation can be evaluated by a sum. 



(B.5) 



{so)=j:T.S{a)P{k,a) 

a k=0 



= Q'i-Q'3 + Q2iQi-Q3)J2Q2 

k=l 

= Ql — Qs ~\ I — Q ■ 
Using the definition of the correlation function (see ([6])) it can be found that. 



(B.6) 



C{r) 



2{Q[ - Q',){Qi - Q3)Q'2Q2-' 
I-Q2 

Q2\Qi-Q3?[ir-l){l-Q2) + Q2]Q 



r-2 



+ g'/(r - i)gr^(i - Q2) + 2Q'2Qi-\i - g'2) + 

1 2 



/2^r-l 



+ 



1-Q2 



Q'2 



(B.7) 
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The final two (constant) terms cancel, and the above formula can be rewritten as, 

f 4Q'2QlQ3 



C{r) 



Q2YQ2 



Q2) + 



/2 



+ 



{i-Q2)Ql' 



{Q2 



which completes the derivation of the correlation function given by f|T3|) . 



(B.8) 



Appendix C. Numerical Support 

The exact expression f|T3|) for correlation function can be supported by means of 
numerical simulations. Since the number of sites in a Cayley tree grows exponentially 
with the number of generations it contains, a useful numerical description of a system is 
computationally expensive on this system. Instead, our numerical calculations have been 
undertaken on a g-regular graph (thin random graph i-e. a graph of vertices, 

in which bonds are placed randomly in such a way that each node has a coordination 
number exactly equal to q. Due to the absence of a boundary, finite size effects are 
much less pronounced in such a system than in a Cayley tree of the same size. This 
graph contains loops, however these loops are small in number, and are irrelevant for 
large enough systems [27j. Therefore, a g- regular graph gives a good approximation to 
a Bethe lattice with the same coordination number. A set of spins is placed on the 
nodes of the g-regular graph. These spins interact according to the Hamiltonian given 
by ([5]) and relax by Glauber dynamics as described in section [2l A similar comparison 
of numerical data for magnetisation versus external field obtained for a g-regular graph 
with analytical results derived for a Bethe lattice has been performed in |8j. 

In our numerical model, each spin in the g-regular graph is assigned a random field 
hi according to a normal probability distribution, hi ~ A/'(0, A^) with zero mean and 
variance A^. All the spins are initally set in the down state, and the system is allowed 
to evolve according to the dynamics described in section [2] until it is in a stable state. In 
order to calculate the correlation function C(r) given by (|6]), a mean (sjSj) is calculated 
for all pairs of spins i and j separated by a chemical distance r. The correlation function 
is then averaged over several (10^) realisations of disorder. 

The r-dependence of the correlation function C(r) presented in Fig. (2) is in good 
agreement with the exact analytical expression derived in this paper. 



Appendix D. Critical Exponents 

The behaviour of the zt-RFIM can be determined by the number of solutions of the 
self-consistent equation G{P*) = (where G{P*) = F{P*) - P*; see (JT])). In the 
presence of external field H, and with random fields distributed according to a pdf 
with zero mean and variance A^, the number of solutions of the self-consistent equation 
depends on the location in parameter space (if. A) [8]. For a Bethe lattice with q = A 
and a normal distribution of random fields, A/'(0, A^), the loci of points corresponding 
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Figure Dl. Locus of points in parameter space [H/ ^1 J) for various number of 
real solutions of the self-consistent equation, G{P*) — 0, with the shaded and white 
regions corresponding to three and one solutions, respectively. At the boundaries 
(continuous and dotted lines), the equation G{P*) = has two roots, one of which is 
a multiple root. The values of {H/J, A/ J) marked by the continuous line correspond 
to the occurrence of the infinite avalanche for increasing external field. These loci 
were obtained for the zt-RFIM defined on a Bethe lattice with g = 4 with random 
fields distributed according to a Normal distribution, A/'(0, A^), in the presence of an 
external field H. The critical point at {Hc/J, AJJ) ~ (1.0, 1.78215895) is marked 
by C. 



to different numbers of real solutions of G{P*) = are presented in figure IDll The 
shaded area (excluding the boundaries) corresponds to the range of parameters where 
the self-consistent equation has three real solutions corresponding to different branches 
of the multivalued function, {m){H,A), for mean magnetisation, (see the part of the 
solid curve in figure |A3] marked by circles). The boundaries correspond to the case of 
two solutions, with the boundary marked by a solid line corresponding to the locus of 
points where an infinite avalanche occurs as external field increases from — oo to oo (see 
arrow in figure R3|) . Outside the shaded region, the self-consistent equation has a single 
solution (see dashed line, dot-dashed line and the part of the solid line not marked by 
circles in figure [A3]). The point at which all three solutions converge is a critical point 
marked by C in figure IDll and represents the boundary between low and high disorder, 
i.e. with and without an infinite avalanche, respectively (see dotted line in figure R3l) . 

In this section, we analyse the behaviour of the correlation length in the zt-RFIM 
around the critical point {A,H) = {Ac, He). In order to do this, we expand (l20l) in 
terms of the differences P* — Pc, H — and A — measuring the distance to the 
critical point, (Pc, Ac, He). The derivative dG/dP* in this equation can be expanded in 
terms of the values P*, H and A near [P^, Ac, He) in the following way. 



dC 

dP* 



P*,HA 



dG 
dP* 



+ 



d^G 
gp*2 



[p* - p:) + 



p:,HcAc 



d^G 
dP*dH 



p;,HcAc 
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+ 



dP*dA 



(A -A, 



dP 



Jf: 3 



*\2 
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where terms of higher order have been neglected because they do not play an essential 
role in the vicinity of the critical point. In the above expansion, we have used the fact 
that the condition for the merging of the three solutions of the self-consistent equation 
is dG/dP* = d'^G/dP*^ = 0. 

The observed critical point is known to be a saddle- node transition [8], meaning 
that the solution, P*, of ([7]) exhibits the known mean-field exponents around this point. 
Explicitly, the behaviour of P* — Pc around criticality with fixed A = Ac can be 
obtained by following the standard method for saddle-point transitions (see e.g. [28]). 
i.e. expanding G{P*) = along the dot-dashed line in figure IDll in terms of P* — Pc 
and H — He- Keeping terms only to lowest order in H — and P* — P^ the equation 
reduces to the following expression, 

1/3 / „ \ -1/3 



P*-Pr = ± 6 



dG 
dH 



d'G 



dP* 



\H-H, 



,1/5 



(D.2) 



PcHcAc 



where ± refers to the sign of H — and the exponent takes its mean-field value, 6 = 3. 
Similarly, expanding G{P*) = along the dashed line in figure IDT] and keeping lowest 
order terms in A — Ac and P* — Pc, the behaviour is found to be, 

-1/2 

1 (- S 1 

PcHcAc 



- 6 



dAdP' 



0, 



PcHcAc 



\A-Ar 



A < Ac 



where the exponent also takes its mean-field value, /3 
and fID.Sp into fID.ip results in. 



(D.3) 

A > Ac 

1/2. The substitution of ([El 



dG 



dP* 



dG 



dP* 



P*,HAc 



P-,HaA/ 



dP*y [dHJ 




\H-Hr 



d^G 
dAdP* 



|A- Ac 



A > Ac 



d^G 



-1 



,(D.4) 
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|A - Ac^^ A < Ac 



dAdP* 

where the critical exponents for the correlation length take the values, u = 1 and 
fi = 2/3. As remarked in section HI these values are different from the mean-field 
exponents, Umf = 1/2 and /imf = 1/3. In order to confirm the values of u and 
yU, we have calculated the spin-spin correlation length for zt-RFIM in the case of 
normally distributed random fields, hi ~ A/'(0, A^) by solving (|8]) and ([7]) numerically 
and using f|TT]) and f|T5|) . The results for the correlation length vs H — He (for fixed 
A = Ac ~ 1.78215895J) and A - Ac (for fixed H = = J) are compared 
with (1D.4P and ( ID.5|) in figure [D2|a) and figure [D2] (b). respectively. A similar check has 
been performed up to g = 10, finding no deviation of the exponents from that for q = 4. 




Figure D2. Spin-spin correlation length for the zt-RFIM in the case of normally 
distributed random fields, hi ^ A/'(0,A^) vs (a) external field, \H — Hc\, at critical 
value of disorder (A — Ac ~ 1.78215895J) and (b) degree of disorder |A — Ad at 
critical value of external field {H = He ^ J). In panel (a), the behaviour on both 
sides of criticaHty is identical (it depends only on \H — Hc\), and is shown by the 
circles. In panel (b), the behaviour is different above and below criticality, and both 
behaviours are shown (squares for A > Ac and circles for A < Ac). The solid line 
in panel (a) and the solid and dashed lines in panel (b) are plots of (|D.4|) and (jD.5|) . 
showing convergence of the exact correlation length to these equations near criticality. 
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